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ABSTRACT 


This study explores the fluxes of heat and salt 
associated with thermohaline staircases in the Beaufort Sea. 
An inverse model is developed to calculate the vertical 
transport of properties in the southern portion of the 
Beaufort Sea directly from observations. The applicability 
of laboratory derived 4/3 flux law is addressed. Three 
formulations of the static advective-diffusive equations are 
discretized on a uniform grid, and inverted using the method 
of Total Least Squares. The first formulation is based on 
the classic flux-gradient model, and is analyzed in both one 
and three-dimensions. The second formulation utilizes 
Turner formulation of the 4/3 flux law, but the formulation 
of the coefficient C (R p ) remains unknown. The third 
formulation includes the full form of the Kelley model, but 
the amplitude of C(Rp) is unknown. Layer averaged heat flux 
through the thermohaline staircases is found to be on the 
order of 1 W/m 2 , which is significantly higher than previous 
studies. This implies that the laboratory formulation of the 
4/3 flux law may require calibration in order to accurately 
represent Arctic conditions. 
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I. INTRODUCTION AND BACKGROUND 

A. ARCTIC CIRCULATION 

Sea ice in the Arctic is melting at an alarming rate, 
much faster than expected (Stroeve et al., 2007). The years 
2007-2009 have seen the lowest summer sea ice coverage in 
the Arctic in the modern era. Ice melt in the Arctic has 
been largely attributed to rising greenhouse gas 
concentrations and the associated warming of the atmosphere. 
The ocean, however, responds to the atmospheric conditions 
on much longer time scales. We are only beginning to 
understand the response of the Arctic Ocean to the 
atmospheric forcing associated with climate change (Stroeve 
& Maslowski, 2008; Dickson et al., 2000). One such response 
is the delivery of significantly elevated guantities of 
warmer water to the Canada Basin and Beaufort Sea. This 
additional warm water, which has taken decades to traverse 
the Canada Basin, adds significant amounts of heat to the 
upper Arctic Ocean, aggravating the ice melt problem. 

The heat stored in the oceans directly influences the 
melting rate of the ice from below. Perovich et al. (2003) 

showed a 1-2 W/m 2 difference between the amount of heat 
necessary to melt the ice from below and the atmospheric 
heat fluxes calculated over open water. They suggest that 
there must be oceanic sources of heat to the ocean-ice 
interface. Understanding the mechanisms that drive the 

heat flux to the ice from the ocean is critical to the 
explanation and prediction of increased sea ice melt. The 
ocean can play a role in melting sea ice via several 
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mechanisms. Ice-ocean interactions include horizontal 
advection of warm water from the Pacific and Atlantic oceans 
under the ice cove; the accumulation of heat due to 
increased solar radiation in ice diminished regions; and 
locally induced upward heat flux into the mixed layer due to 
upwelling, topographically controlled flow, and eddies 
(Maslowski et al. , 2009). 

This thesis will focus on the role of double diffusion 
in oceanic heat transfer, which is a subject of growing 
interest in oceanography and Arctic sciences. It is argued 
that double-diffusive convection may also play a significant 
role in transferring heat upward into the Arctic mixed 
layer. The heat transported toward the ocean surface via 
diffusive convection may be a critical component to 
understanding the connection between ocean dynamics and sea 
ice extent in the Arctic. 

B. DOUBLE-DIFFUSIVE CONVECTION 

Double diffusion is defined as a set of processes 
related to the difference in molecular diffusivities of heat 
and salt. It is known to influence vertical mixing in the 
ocean, and comes in two forms: salt fingering and diffusive 
convection. Salt fingering is prevalent in the mid¬ 
latitudes, while diffusive convection is the form of double 
diffusion seen in the Arctic. 

The term diffusive convection is used to describe two 
types of motion that arise in fluids with a negative 
vertical gradient of temperature and salinity. The first is 
characterized by oscillations in smooth temperature-salinity 
gradients. The other, diffusive layering, operates in a 
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thin interface, which separates two homogeneous layers. Heat 
and salt are transferred across the diffusive interface. 


The oscillatory regime of double diffusive convection 
is described in Figure 1. 


cold, fresh 



Figure 1. The oscillatory regime of double diffusive 

convection 

Imagine a parcel of water in a lower layer that is 
displaced upward. This parcel is warmer and more saline 
than the layer of water above. As the parcel enters the 
layer above, it immediately begins to lose heat, via 
molecular diffusion, to the surrounding upper layer at a 
rate nearly 100 times as fast as it loses salt. The result 
is cooling of the parcel, while it maintains its salinity. 
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The parcel then begins to descend, due to its increasing 
density, as it cools. The denser parcel begins to gain heat 
in the lower, warmer layer, again via molecular diffusion, 
and once it gains a sufficient amount of heat its density 
decreases and the parcel begins to rise again. This mixing 
results in the formation of layers of near constant 
temperature and salinity. It is these layers that define 
the second form of diffusive convection, diffusive layering. 
This process is known to be prevalent in the Beaufort Sea. 



Potential Temperature fC) 



Figure 2. Temperature - Salinity (From Timmermans, 2008) 


Diffusive layering occurs in a 
where, in the case of the Arctic, cold 
warm and salty water, as seen in Figure 


stable environment, 
fresh water overlies 

2 . 
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Figure 3. Diffusive Layering - Fluxes of heat and salt are 
transported upward across the thin layer, causing 
convection and mixing within the homogeneous layers. 


Diffusive layering is manifested in the Arctic 
thermocline by a series of homogeneous layers of near 
constant temperature and salinity separated by a thin 
diffusive interface. This natural occurrence is often 
referred to as a thermohaline staircase. The physical 
processes associated with the maintenance of this structure 
(Figure 3) can be described as follows: Heat and salt are 
transferred upward through the diffusive interface. The 
water just above the interface gains sufficient heat to 
begin upward motion due to excess of buoyancy. This water 
begins to rise within the layer, driving the convection that 
ultimately maintains the homogeneous properties of the mixed 
layer. 
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The data shown in Figure 2 are representative of the 
majority of the data collected by the Woods Hole 
Oceanographic Institute (WHOI) Ice Tethered Profiler (ITP) 
program. The warm water seen from ~300 meters and below 


originates in 

the 

Pacific 

and 

Atlantic 

oceans 

and is 

advected into 

the 

Beaufort 

Sea 

via the 

Arctic 

current 

system. The 

cold 

water above 

is a result of 

downward 


temperature flux associated with the cold surface and air- 
sea-ice interaction. The less saline, or relatively fresh, 
surface water is a result of annual ice melt and fresh water 
flux associated with river runoff. These conditions, and 
the lack of vertical mixing below the surface mixed layer, 
provide an ideal environment for the formation of double 
diffusive staircases. 


C. EARLIER ESTIMATES OF HEAT FLUX 


The fluxes associated with double diffusion have been 
studied extensively. It is generally accepted that in 
systems characterized by active diffusive layering, the 
interaction between layers is controlled by the flux of salt 
and heat across the diffusive interfaces. Based on a set 
laboratory experiments the balance of temperature variation 
between adjacent layers and heat flux (Turner, 1973) was 
proposed as follows. 


4 

aF T = Aj -[alsiy (1) 


, _ /?A 5 
' p a AT 


( 2 ) 


This relationship (1) has since been known as the 4/3 
flux law. Ax has the dimensions of velocity and is a 
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function of the density ratio (R p ) . F T is temperature flux, 
a is the thermal expansion coefficient, and (3 is the saline 
contraction coefficient. 

The 4/3 flux law states that the fluxes associated with 
double diffusion are controlled entirely by the variation in 
temperature or salinity across the diffusive interface and 
is independent of layer thickness. The lack of sensitivity 
of fluxes to vertical scales is critically important to our 
observation-based study because the thickness of an 
interface is often much smaller (<25cm) than the resolution 
of the sensor, whereas, a difference in mean temperatures of 
each layer can easily be measured. Dimensional analysis 
suggests that the 4/3 flux law can be written as follows: 

i 

aF T =c(R p )\^- ) \ \a\Tf (3) 

V v J 

C (R p ) is an unknown constant, g = 9.8 m/s, v = 1.8X10 -6 
m 2 /s is the kinematic viscosity, k = 1.4X10 -7 m 2 /s is the 

molecular diffusivity of heat. 

Marmorino and Caldwell (1976) examined the 4/3 flux law 
in depth. They discovered an empirical formulation of the 
coefficient C(R p ) using the following equations: 

if = 3.8/y (4) 

H s P 

»„= (0.085)*, 

-^- = 0.101exp|4.6exp^-0.54(/? /i -l) j (6) 

H sp 
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C(/?J = 0X)0859exp(4.6exp[-0.54{/^-l}]) 


(7) 


Using Huppert's formula (4), they were able to 
determine a relationship between the measured heat flux (H) 
and (5), the theoretical heat flow through a non-deformable 
interface (H sp ) . They proposed a new formulation also based 
on laboratory experiments (6), which ultimately yielded a 
value for C (R p ) . 


The coefficient C(R p ) has since been scrutinized, most 
notably by Kelley (1990) . In his 1990 paper, Kelley 
proposed a new value for C (R p ) based on what he referred to 
as a full collection of laboratory measurements. He 
analyzed laboratory experiments for l<ih,<10, while Turner 
(1973) only looked at values of R p up to 7. 


C(R P ) = 0.0032 exp 


4.8 

0.72 




) 


( 8 ) 


Equation (8) is the result of Kelley's laboratory 
experiments, and along with (7), the two have been used 
interchangeably as a standard method for calculating fluxes 
using the 4/3 flux law. However, Kelley suggested that 
perhaps the 4/3 exponent is incorrect and should be changed 
to 5/4. Additional experimentation is required to support 
this proposition. 

Padman and Dillon (1987) showed that using the 
Marmorino and Caldwell formulation of the 4/3 flux law, they 
could calculate the fluxes in the Beaufort Sea. Their 
results showed that those fluxes were 0.02 — 0.1 Wm -2 . A 
more recent study (Timmermans, 2008), also based on 
extrapolation of laboratory derived results, estimated the 
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heat flux associated with the diffusive staircases as 
measured by the Ice Tethered Profilers to be 0.22 + /- 0.10 
W/m 2 . 

Wilson (2007) compared the two formulations (7) and (8) 
using data collected via the Ice Tethered Profilers, a set 
of moored CTDs in the Beaufort Sea, which will be discussed 
in detail in Chapter II. She found that in order to fit the 
data collected in the Beaufort Sea to these formulae each 
needed to be multiplied by a transfer coefficient. She 
suggested that the laboratory results of both Kelley (1990) 
and Marmorino and Caldwell (1976) cannot be applied to 
actual data collected in the field, and perhaps the 
laboratory measurements require calibration prior to 
application in the Arctic. She concluded that the average 
heat flux in the Beaufort Gyre in excess of 1 Wm~“. This is 
significantly larger than estimates using the original 4/3 
flux law. 

Caro (2009) used a set of high-resolution numerical 
experiments to validate some of the conclusions made by 
Wilson. He used two-dimensional numerical models to 
estimate the heat fluxes in the Beaufort. His efforts 
resulted in fluxes on the order of 1 W/m". 

Timmermans' findings differ from Wilson's (2007) and 
Caro's (2009) by an order of magnitude, and in some cases 
Padman and Dillon's differ by two orders of magnitude. Such 
a wide range of estimates in a region critically important 
for climate variability provides the impetus for this study. 
This thesis attempts to resolve the controversy in 
assessment of the vertical heat transport by applying 
inverse modeling techniques to observations. 
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D. INVERSE MODELING 

The study of the oceans traditionally follows two 
distinct directions. One involves numerical modeling and 
theoretical manipulation of the equations governing dynamics 
and thermodynamics of the sea. The other is observationally 
oriented, where real measurements are discussed in the 
context of qualitative physical principles. 

As oceanography matures into a quantitative science, it 
becomes increasingly desirable to connect the observational 
and modeling aspects of the field. Laboratory experiments 
often attempt to bridge this gap between hard facts and 
theoretical deductions. Laboratory researchers strive to 
replicate the ocean environment as closely as possible, but 
exact correspondence is seldom achieved, given the 
substantial differences in the conditions and spatial scales 
between the laboratory and nature. 

Another more recent approach, which makes it possible 
to connect data with theory, involves inverse modeling. 
Inverse models attempt to directly calculate model 
parameters from observations by minimizing errors associated 
with observational measurements and fitting those 
measurements to the governing equations. In our case, 
temperature, salinity and depth measurements were taken by 
the Ice Tethered Profilers, and velocity and diffusion 
coefficients are calculated from those observations using 
the method of Total Least Squares. 

v.(vr)+dl = L( Fr ) (9) 

OZ OZ 

V«(V5 ) + ^ = A( F ) (10) 

v 7 dz dz y J 
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( 11 ) 


V.V + —= 0 

dz 

The inverse model in this study is based on the static 
advective diffusive (9), (10) and mass continuity equations 

(11) in divergence form. F T , and F s is the temperature or 
salinity flux. The formulation of F T and F s is central to 

this study. In this work we examine flux-gradient 

formulations in which F T , and F s are assumed to be 
controlled by the background temperature and salinity 

gradients, as well as the interfacial flux laws. Results of 
both are qualitatively consistent, but differ in details. 


E. 


UNRESOLVED ISSUES 


This thesis attempts to quantify the heat flux 
associated with diffusive convection independent of the 
extrapolation of laboratory derived laws. This information 
is used to address the following questions. 

1. What are the typical vertical heat and salt fluxes 

in the diffusive staircases in the Beaufort Gyre? Vertical 
heat fluxes of 0.02 — 6 Wm -11 represent an extremely wide 
range making it difficult to assess the large-scale 
consequences of heat transfer in the Arctic. If the value 
is less than 1 Wm” 2 , then it is likely that the heat 

delivered upward vis-a-vis diffusive convection is not a 
major contributor to sea ice melt. However, if Wilson 
(2007) is correct, then diffusive convection is sufficient 
to affect the pattern on stratification and ocean climate in 
the Arctic. 

2. Can the heat flux through the diffusive staircases 
significantly contribute to the heat imbalance associated 
with the melting of the polar ice cap? 
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3. What are the relative roles of the vertical and 
lateral transport of heat in Beaufort gyre? Is it possible 
to capture the zero-order dynamics at play by a one¬ 
dimensional horizontally averaged model? 

4. Is it possible to accurately evaluate fluxes 
associated with double diffusion by inverse modeling 
techniques? Lee and Veronis (1991) showed that it is 
possible to quantify the fluxes via an inverse model using 
the method of Total Least Squares in the tropical Atlantic 
Ocean. A similar inverse model of the advective-diffusive 
equations is applied to data collected from 2004-2009 via 
the Ice Tethered Profilers operated by the Woods Hole 
Oceanographic Institute. Using methods described in Lee and 
Veronis (1991) the inverse technique is implemented to 
determine eddy diffusivity coefficients using the method of 
total least squares. A comparison to the estimates using 
the 4/3 flux laws described in Timmermans (2008) and Padman 
and Dillon (1987) is discussed along with a comparison to 
Wilson (2007), additionally the results of this study area 
compared to the model results from Caro (2009). 

This thesis is organized as follows: Chapter II reviews 
the data from the Ice Tethered Profilers and some of the 
unique processing required for obtaining the appropriate 
data set for the inverse calculation. Chapter III offers a 
detailed description of the mathematical foundation of the 
inverse model as well as its application to the ITP data 
set. Chapter IV presents the results of four models, the 
one dimensional (vertical only) model, which uses an average 
single profile to determine fluxes, and three variations of 
a full three dimensional model. A discussion of the 
findings and conclusions is in Chapter V. 
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II. DATA AND DATA PROCESSING 


A. ICE-TETHERED PROFILER DATA 

The Woods Hole Oceanographic Institute (WHOI) Ice 
Tethered Profiler (ITP) program was established in 2004 to 
provide timely and accurate oceanographic measurements in 
the Arctic Ocean. Historically, access to the ice covered 
Arctic Ocean has been limited. The ITPs assist in assuring 
consistent measuring capability in a data-sparse 
environment. Between 2004 and 2006, the first six ITPs were 
deployed in the Beaufort Sea. Data from ITPs 1-6 are used 
in this study. As of July 2009, these six sensors collected 
over 7000 profiles, all of which have been made available to 
the public. 

The ITP is a moored conductivity, temperature, and 
depth (CTD) sensor. Each one is designed to move with the 
ice it is moored to and has the potential lifespan of 
approximately three years. The ITP consists of a small 
platform on top of which rests a buoy. A 10-inch hole is 
drilled through the sea ice and the instrument is inserted 
in the water from above. The instrumentation is tethered to 
the surface float via a wire rope that is weighted on the 
bottom. This rope is approximately 800 meters in length, 
and the instrument traverses the wire rope several times 
daily. 
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Figure 5. Location of ITP 1-6 

Figure 4 illustrates the ITP system and Figure 5 
presents the trajectories of the first six profilers 
followed over their lifespans. The data used in this study 
are labeled "Level 3 Archive Data" by WHOI. This form of 
data represents the best possible representation of the 
ocean properties as measured by the ITPs. Corrections have 
been applied that account for sensor response, regional 
conductivity variation, and quality assurance. Each data 
file is formatted to be compatible with MATLAB. Each 
profile was retrieved via 
ftp://ftp.whoi.edu/whoinet/itpdata/ and processed further as 
described below. 
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B. DATA PROCESSING 

The corrected ITP data described above was arranged by 
profiler number 1-6 and processed together as one data set. 
Processing began with the corrected ITP data obtained from 
the WHOI ITP Web site. Wilson (2007) suggested that 
anomalous results were obtained from the profiles 
corresponding to the downward data collection for each cast. 
This is possibly due to sensor contamination associated with 
the main ITP sensor platform below the sensor suite. As the 
profiler is lowered into the water the platform interacts 
with the water just below the sensors, affecting the 
measurements. Based on this suggestion, all downward 
collected profiles were removed from the total dataset, 
effectively reducing the amount of data by half. 

Diffusive convection occurs in the thermocline region 
of the Beaufort Sea between approximately 200 meters and 400 
meters depth, with the majority of the diffusive layers in 
the upper 100 meters of the thermocline. The data were then 
reduced to only those data from 200 db pressure to 300 db 
pressures, which correspond to roughly 195 meters to 295 
meters depth. All temperatures were converted to potential 
temperature. Table 1 summarizes the remaining data. Each 
profile consists of temperature, salinity, pressure, and 
depth measurements. 
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UP 

Max # of Obs in 
each Profile 

# of Profiles 

1 

1051 

1022 

2 

520 

122 

3 

1987 

766 

4 

1917 

349 

5 

2327 

547 

6 

1710 

667 


Table 1. Number of Observations (individual data points) 
for Each Profile. Each profile has up to the maximum 
number of observations indicated. Not all profiles 
have the maximum number of observations. 


These data were further reduced by determining the 
values of temperature and salinity at the center of each 
layer. This was accomplished by searching the data for the 
local maxima of vertical temperature gradient. This local 
maximum corresponds to the interface between each layer. 
Using the position of each interface, the center of each 
layer is determined to be exactly half way between two 
adjacent interfaces. The temperature, salinity and depth of 
these data points is recorded and marked. Table 2 shows the 
number of individual diffusive layers found in the data 
organized by profiler. 


UP 

# of Layers 

# of Profiles 

1 

48 

1022 

2 

39 

122 

3 

51 

766 

4 

54 

349 

5 

51 

547 

6 

62 

667 


Table 2. Number of Layers for Each Profile. 
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Figure 6. Temperature - Salinity plot for ITPs 1-6. 

The diffusive layers next must be classified 
independently of the particular cast. Timmermans (2008) 
showed that each diffusive layer can be identified via a 
plot of potential temperature versus salinity. Figure 6 is 
a Temperature-Salinity plot for each ITP used in this study. 
The vertical clusters at nearly constant salinity are 
associated with a horizontally homogeneous layer, a 
diffusive layer. Due to the remarkable lateral coherence of 
layers in the thermohaline staircase, it is possible to 
consolidate the data into one dataset by identifying 
pronounced layers that appear in all six datasets. 
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Figure 7. Histogram of data: Number of data points for a 

given salinity value. 

Using a moving window across salinity space that 
corresponds to a 0.001 PSU increment from 34.3 PSU to 34.6 
PSU we estimated the concentration of data points as a 
function of salinity, as shown in Figure 7. The blue curve 
is the original histogram. The red curve is a 5-point 
moving average, which is used to identify the individual 
layers. Noting that the concentration of points is highest 
at the center of the layers and lowest at the interface, we 
assume that the troughs in the histogram represent the 
boundary (interface) between each layer. The black dashed 
lines represent the inferred interface locations. 
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UP 

# of Layers 

# of Profiles 

1 

19 

1022 

2 

19 

122 

3 

19 

766 

4 

19 

349 

5 

19 

547 

6 

19 

667 

ALL 

19 

3452 


Table 3. Total number of layers and profiles in final data 

set. 

These values are used to enumerate layers in all 
profiles. The result is a total of 19 layers common to each 
profiler (Table 3) . The data between adjacent troughs 
(dashed lines in Figure 7) are considered to be in the same 
layer. Thus, layer 1 contains data only from the region 
between the first two troughs on the left of the graph. 
Layer 1 corresponds to the shallowest layer in the 
thermocline, and layer 19 corresponds to the deepest (as 
salinity increases downward). It is important to note that 
diffusive layers in the Beaufort Sea are not limited to 
these 19. These are the layers that are common to all 6 
ITPs used in this study. 

The highlighted portion of Figure 7 shows the 
difficulty in capturing each layer accurately. The blue 
line clearly shows a highly variable structure in this 
region, and our method for determining the interfaces 
between adjacent layers is not precise enough to account for 
this variation. As a result, we observed some anomalous 
results in the inverse calculations. 
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Figure 8. Locations of each profile. 

Ultimately, each profile yielded a specific value of 
temperature, potential temperature, salinity, layer 
thickness, depth, pressure, and latitude and longitude for 
each of the 19 layers. Figure 8 is a plot of the locations 
of all 3452 profiles. The next step involves interpolation 
of the data onto a structured grid for use with the inverse 
model. 

C. INTERPOLATING DATA TO A GRID 

The inverse model used in this study requires the data 
to be applied to an equidistantly spaced grid in both 
horizontal directions. This is accomplished using the 
MATLAB routine "griddata," which linearly interpolates the 
data to a pre-determined set of grid coordinates. The 
result is a dataset that is compatible with the inverse 
model. The horizontal length scale is 15 km, both zonal and 
meridional, unless otherwise noted. The vertical length 

scale is defined by the layer thickness at each location. 
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3D Model Domain 



Figure 9. Location of the three-dimensional dataset. 

The grid used for the three dimensional model has the 
dimensions 7X7X19. Figure 9 shows the location of the grid. 
This location in the southern Beaufort Sea was chosen 
because of the high density of original profile locations, 
as well as the relatively even distribution of profiles 
throughout the grid domain. 
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Figure 10. (a) Salinity surfaces for 19 layers used in the 

inverse calculation, (b) The corresponding potential 

Temperature surfaces. 

The structured data are shown in Figure 10, which 
reveals nearly constant temperature and salinity values in 
each of the diffusive layers, and the lateral coherence of 
those layers over of the model domain. For both Figures 10a 
and 10b, the layers are arranged from the deepest to the 
shallowest from top to bottom. The salinity values for each 
layer are nearly constant spatially. The temperature values 
vary slightly, but the homogeneity is evident. 
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III. INVERSE MODEL 


A. MODEL FRAMEWORK 


In order to apply (9), (10) and (11) to measurements 
they must be discretized in space. This is achieved by 
integrating vertically from the bottom (b) to the top (t) of 
each layer, assuming vertically uniform temperature and 
salinity within each layer, consistent with characteristic 
properties of diffusive convection in thermohaline 
staircases. 


V* 


V* 


f— \ 
XT 
V J 
f~ \ 

xs 


+ [wT-F T l=0 
+ l"S-F 5 1=0 


V J 

vnv+[w]' =o 


( 12 ) 

(13) 

(14) 


V represents the vertically integrated horizontal 
velocity along the layer, and w represents the vertical 
velocity across the interface. The vertical influences into 
and out of each layer are described by the diffusive and 
convective terms in (12) and (13) (Lee et al. , 1991) . 


The inverse model is based on the optimal solution to 
the overdetermined linear system of equations Ax=b. 
Literature refers to A as the data matrix and b as the 
observation vector. We will follow this same convention in 
this study. The solution is x, a vector with the same 
dimensions as b where each component of x represents the 
coordinates of b in terms of the columns of A. 
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Figure 11. The seven-point grid network for finite 
differencing. (After Lee & Veronis, 1991) 


To formulate the problem in terms of the algebraic, 
rather than the differential system the governing equations 
are discretized in three dimensions using centered finite- 
differencing on a seven-point grid (Figure 11) . Each cell 
in the grid has the measured tracer value, temperature and 
salinity, in the center, u and v on the lateral boundaries 
and w, F t and F s on the vertical boundaries. 
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The model discretization is employed exactly as 
described in Lee and Veronis (1991). Equations (15), (16) 
and (17) are the discretized forms of the advective 
diffusive equations and the continuity equation 
respectively. The indices i, j, k represent the location of 
the tracer for that qrid cell, where i increases eastward, j 
increases northward, and k increases downward with depth. 
The ^ values correspond to the boundaries of the grid cell. 
Velocities and diffusive fluxes into and out of each cell 
are calculated at the boundary between grid cells. Thus, 
the coefficients are weighted averages of thickness and 
tracer values at each corresponding boundary. Slight 
variations are used in the three-dimensional models as we 
explore different formulations of flux calculations. These 
variations are described in Chapter IV. 

Notice that the integrated forms of the original 
equations are homogeneous. Thus, the model will have an 
exact trivial (zero) solution. To avoid this unphysical 
solution one must introduce a known, or estimated, quantity 
that can be used to populate a non-zero observation vector 
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b. This known quantity will be referred to as the reference 
unknown x r . To illustrate the concept of a reference 
variable consider a homogeneous system (b=0) of four 
equations with three unknowns. 

1 + £ 2 — £ 2 + S 

A-s 2 + s 1 — s 
8 + £ \-s 5 + £ 

6-s 5 + £ 3 + s 

Let us assume that the value of xi is known or can be 
readily estimated with some degree of certainty. All of the 
equations can be divided by xi=x r . This effectively forces 
an inhomogeneous term where bAO and creates a non-trivial 
solution to the problem. 

2 -a 3 + £^| f 1 + £ 

2 + s 1-s (x 2 lx^\__ 4-f 

\-£ 5 + £ lv X 3 / X l J 8 + ^ 

5 + ^ 3 + sj l6-f 

Lee (1991) showed that this is an effective method to 
produce the inhomogeneous terms that are necessary for the 
inverse calculations. The model, however, is very sensitive 
to the choice of reference variable. This sensitivity is 
explored in Chapter IV. 

Another essential step to obtain a reliable solution to 
the problem is to ensure that the system of equations is 
overdetermined, i.e., that we have more equations than 
unknowns. This acts to smooth out the errors in the data. 
Each grid cell contains three equations, two tracer 
(temperature and salinity) equations, and continuity. As 
written, there are three equations and five unknowns, u, v, 
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w, F s , and F T . Two important assumptions are made to ensure 
that there are fewer unknowns than equations. The first 
assumption is that w, F s and F T are constant at each 
interface, and the second is that w at the top layer is 
equal to Ekman pumpinq velocity (w e 5 x 1CT 7 m/s) at the 
top of the thermocline (Yang, 2006). This ensures that the 
system is overdetermined. In our case, we have 1275 
equations and 1073 unknowns 

Finally, the data matrix A and the observation vector b 
are built using the coefficients in (15), (16) and (17). 

This is accomplished by populating a matrix with the 
coefficients of the model unknowns in the correct location. 
Each equation is represented by a specific row and each 
unknown by a specific column. The result is a very sparse 
matrix, where most of the values in the matrix are zero. 



coeff (mj ) 

coeff (u 2 

) o 


0 -A 



0 

coeff (u 2 

) coeff (u^ 

) 0 

0 


A = 

0 

0 

coeff' ( u 3 

) coeff (u 4 ) 

0 

(20) 


0 

0 

0 


0 



, 0 

0 

0 

0 coeff (u n ) 

coeff (u n+ i)--y 



Equation (20 

) is an 

example of the 

structure 

of the 3- 

dimensional 

form 

of the 

data matrix A, 

and n corresponds to 


the number of u unknowns. A similar pattern continues for 
the coefficients of v, w, F T , and F s completes the matrix. 
In our case A has the dimensions 1275 rows (equations) X 
1073 columns (unknowns). 

The reference unknown and its corresponding 
coefficients are divided through as shown in (19), and the 
vector b is obtained by subtracting the coefficients for the 
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reference variable from both sides of the equation. The 
basic framework for the TLS method and the inverse 
calculation is now complete. 

B. THE "BEST" SOLUTION TO THE OVERDETERMINED SYSTEM 


The methods 

A e R"' xn and the 
relatively well 
therefore Ax=b is 
has more equations 
exact solution, 
formulate a "close" 

Least squares fit 



a i 


data matrix 

given, are 
m > n and 
e., one that 
is generally no 
one attempts to 
solvable. 


Total least squares fit 



“i 


for solving Ax=b, where the 

observation vector be R n are 
known. Specifically, when 
an over determined system; i. 
than unknowns, there 
In such cases then, 
problem that is exactly 


Figure 12. Least squares and total least squares fits of a 
set of m=20 data points in the plane and just one 
unknown x. ai and bi represent the components of an 
example data matrix a and an observation vector b. In 
the least squares fit the error only exists in b and in 
the total least squares fit errors exist in both a and 
b. (From Markovsky et ai., 2007) 
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Least Squares Fit 


Total Least Squares Fit 




Figure 13. Least Squares and Total Least Squares fits of a 
plane A and vector b. The plane A is a representation 
of all possible values for Ax. Ultimately only one 
vector Ax represents the closest approximation of b. 
(a) LS minimizes the vector -r, where -r only 
represents error in the vector b, by projecting b onto 
the plane A. (b) TLS minimizes -[E|r] which represents 
the error in both A and b, needed to project b onto the 
perturbed plane A. Both the data matrix and the 
observation vector are corrected to find an exact 

solution. 


The classical method of Least Squares attempts to find 
a "best" solution by assuming that no exact solution exists 
entirely due to error in the observations. A in theory is 
assumed to be precisely correct and some small correction r 
exists such that Ax=b+r is exactly solvable. (In this case 
-r represents the errors in the observation vector.) The 

solution of the problem is a vector xeR" corresponding to 
the minimum correction required to make the problem 
solvable, i.e. min||Ax-b|| = minllr II . Geometrically, the 

xe R" 11 112 xe«" 11 112 

solution represents either minimizing the sum of the squared 
vertical distances from the data points to the fitting curve 
(Figure 12a), or the projecting the data vector onto the 
plane spanned by the columns of A (Figure 13a). 

The total least squares method is an extension of the 

classical Least Squares technique that takes into account 
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errors in both the data matrix A and the observation vector 
b (Golub et al., 1980) . The method assumes that only 
(A+E)x=b+r is exactly solvable. E represents the correction 
in the data matrix A and r represents the correction in the 
observation vector b. 

I m n 2 

|[E|rhvE2>lh,| (21) 

V «'=1 j =1 

The solution corresponds to making the minimum 
correction in the sense that min ||[E I rill , where IHI denotes 

b+reR< A+E > 11 ' llF 11 llF 

the Frobenius norm, and E e R" Lxn and r e R" . The Frobenius 
norm is defined as square root of the sum of the absolute 
squares of the elements of the matrix (21). 

Geometrically, the total least squares solution can be 
described as either minimizing the sum of the squared 
orthogonal distances from the data points to the fitting 
curve (Figure 12b), or the projecting the data vector onto 
the plane spanned by the columns of A+E (Figure 13b). 

The Total Least Squares and Least Squares methods are 
shown schematically in Figure 12 and Figure 13. The 
difference in the solutions for the two methods amounts to 
changing both the slope of the line and the intercept 
(Figure 12), or adjusting both the plane A and vector b 
(Figure 13) to find a solution, represented by the fit line 
in Figure 12 and the colored vectors on the plane in Figure 
13. 

There are two other important characteristics of the 
TLS solution to consider. First, the method implicitly 
assumes that the values of the data matrix A and observation 
vector b are roughly the same. In Figure 12, the range of 
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and b± are roughly the same, approximately -2 to 2. This 
is important to limiting the sensitivity of the solution 
since the method of TLS makes absolute, not relative, 
changes to A and b. If the range of each input matrix is 
nearly the same, then the norm of the solution vector x will 
be approximately 1. If there are one or more orders of 
magnitude difference in the range of the data matrix A and 
observation vector b then it is likely that the solution to 
the total least squares problem will be greatly different 
from the least squares problem. A large difference between 
TLS and LS can be interpreted as a warning sign that the 
solution may be physically inconsistent, as illustrated in 
Figure 14. A difference by a factor of four in the norms of 
A and b can result in dramatic differences in solutions. 


.east Squares Solution 



Total Least Squares Solution on the Ongiral data 



oast Squares Snlirtinn 



Total Least Squares Solution on the Ongiral data 



Figure 14. The importance of similar norms for A and b. (a) 
A properly scaled problem, where A and b have the same 
norm. (b) The norms of A and b differ by a factor of 

four. 


The second point is that the method implicitly assumes 
that the relative errors associated with the data matrix A 
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and observation vector b are comparable. A large error in 
either input matrix will skew the TLS solution relative to 
the pattern realized in nature. 

An accurate and reliable estimation of E and r are 
paramount to the TLS solution. The value of the correction, 
[E|r] , tells us how far to move or perturb the space of 
[A|b] to obtain the solution x. This is accomplished via 
the Singular Value Decomposition (SVD). 


This process begins with the augmentation of the data 
matrix A with the observation vector b giving a new matrix 
C=[A | b] . The matrix C is factored into three distinct 
components orthonormal matrices U and V and a diagonal 
matrix E. 


C = [A lb] = UXV T 


(u (1> u ,2) ••• u (n+1) ) 


cr, 0 • • • 0 

0 <J 2 • • • 0 


v° 0 ••• ° n J 


f v d)T A 

„(2)T 


,(n+l)T 


: V (1)T + a 2 U ( 2 ) V (2)T + • • • + a„U (n+ 1 ) V (n+1)T 


( 22 ) 


(23) 


(24) 


Eguations (22), (23), and (24) are three formulations 
of the SVD. The result is a series (24) where the first 
term represents the largest singular value of the data, an 
analogue to the largest amount of energy, and the last term 
represents the least amount of energy in the system. This 
is very similar to Fourier series where the terms are 
multiple harmonics from lower freguencies to higher 
frequencies. SVD simply goes from higher energy to lower 
energy. 
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An important distinction of SVD is related to the 

ordering of the coefficients oi. The magnitude of each 
value of Oi decreases with i, and therefore the most energy 
in the system is contained in the first component and the 

least energy is contained in the last component of the 

system. 

As noted before the TLS solution corresponds to the 

correction min ||[E|rl|| where ||*|| denotes the Frobenius 

b+rER (A+E > 11 llF 11 llF 

norm. Given the SVD, the Frobenius norm can be easily 
calculated. 

II*L = \l a i +a 2 + '" + a n (25) 

The Frobenius norm (25) differs greatly from the 2-norm 
used in the Ordinary Least Squares calculation. Recall, the 
solution to the Least Squares problem, minIIAx-bll = minllrll , 

xe R" 11 112 xeR" 11 112 

where, ||C|| = o l . The 2-norm only considers the singular 

value with the highest energy (oi), whereas, the Frobenius 
norm considers a range of singular values. Thus, the 
Frobenius norm contains a more inclusive representation of 
the system. 


C = CS-jiT v + o- 2 u (Z V ZM H— + cr„u n) v 


(n).,(n)T 


A = -a „, u ( " +1) v (n+1)7 ’ = -Cv' 


(n+l) y (n+l)T 


cr n+l = mm 

rank(C+A)<n+l 


(26) 

(27) 

(28) 


Given (24) , the closest (in a |»| sense) rank n C to 

C is given by (26), which is simply (24) with the last term 
dropped. Therefore (27) and (28) are true. 
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X 


(29) 



Consider C= [A | b] , which has rank (n+1) . A is rank n 
and when augmented with b the rank becomes (n+1) . The TLS 
problem can be reformulated as (29), which implies that 

C = [A + E |b + r] is no more than rank n, because C is equal to 

C with the last term removed. The last term in C is used to 
formulate the correction matrix A. 


Ultimately, what is sought is a matrix A = [E | r] e R iraU,+l) 
that describes the minimum perturbation to the data matrix 

A, and the observation vector b such that C is rank n and 
the solution for (29) exists. 

It follows that the solution can be determined via the 


last (n+1) 


column of V, defined as 



where yeR", 


and a^O. If x = - [y/a] , and A = [E|r] = -Cw T then (29) is 
solved exactly. 


C. SCALING THE MATRIX 


Ensuring that each equation is treated with equal value 
in the model requires that we normalize each row. This has 
two purposes: one is to prevent the appearance of spurious 
solutions, and the other is to up weight the more reliable 
equations (Lee, 1991). In this study, we divide each row by 
its Frobenius norm (for vectors the Frobenius norm is 
equivalent to the 2-norm), so that all of the equations have 
the same norm in the system. 
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It is sensible to presume that the continuity equation 
is less prone to error than either tracer equation, since 
there are no tracer values in continuity. For this reason, 
we multiplied the continuity equation by 10, effectively 
statinq that continuity is more reliable than the other two. 
The value of the multiplier is important, but choosing 10 or 
100 ultimately does not change the outcome. We simply give 
continuity a higher weight (relative importance to the 
problem) than the tracer equations. 

D. PARAMETERIZATIONS OF THE VERTICAL FLUXES 

Four models are discussed in this work, each having a 
different formulation of the fluxes (F s , F T ) . 

F T = K T — (30) 

dz 

aF T =c(R p )(aM) m (31) 
olF t =BOO(R p )(aAT) 413 (32) 

Equation (30) is the classical form of flux that is 
dependent on the gradient of the tracer (temperature or 
salinity). This form is explored in both a one-dimensional 
model and a three-dimensional model. 

Equations (31) and (32) are based on Turner's 4/3 flux 
law and invokes a more specific formulation proposed by 
Kelley (1990). Equation (31) is explored where the 
formulation of C (R p ) is unknown, while (32) has the full 
form of C(R p ) and the amplitude is unknown. 

In Chapter IV, we present the one-dimensional model, 
which corresponds to the (30) formulation of the fluxes, and 
the three three-dimensional models: Model 1, Model 2, and 
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Model 3, which correspond to the flux formulations in 
equations (30), (31), and (32), respectively. 

Recall, the solution vector x&R" is composed of n 
elements that represent the unknowns (u, v, w, F s , F T ) in 
the discretized model. Each element in x must be multiplied 
by x r to recover the correct value of the unknown. The 
output of the model is the vector x. 
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IV. MODEL RESULTS 


The previous chapters described the framework necessary 
to implement the model. In this chapter, the results of the 
four different models are presented. As discussed in 
Chapter III, some of the elements of the model framework are 
adjusted to better understand the physical processes at work 
in the diffusive staircases. Each model will be described 
independently. 

A. ONE-DIMENSIONAL MODEL 

The one-dimensional model was developed as a precursor 
to the larger, more dynamically inclusive three-dimensional 
model. The one-dimensional model allowed us to understand 
how the model worked and what type of output to expect. It 
also allowed us to observe the role of horizontal advection 
in the flux calculations. 
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(33) 


(34) 


(35) 
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In this case, we neglect all horizontal velocity 
components and consider a system where vertical advection of 
salt and heat are balanced by diffusion. Equations (33), 
(34) and (35) are used for the inverse one-dimensional 
calculation. Depth is measured from the sea surface and 
therefore vertical velocity is positive downward. 


Layer 

Potential 

Temperature 

(degC) 

Salinity 

(PSU) 

Layer 

Thickness 

(m) 

Layer 

Depth 

1 

- 0.6533 

34.3128 

1.7605 

220.5104 

2 

- 0.6065 

34.3307 

1.9772 

222.3872 

3 

- 0.5589 

34.3484 

1.9700 

224.4895 

4 

- 0.5080 

34.3654 

2.0280 

226.6684 

5 

- 0.4513 

34.3847 

2.1539 

229.5241 

6 

- 0.4074 

34.4006 

2.1856 

232.3383 

7 

- 0.3687 

34.4140 

2.2638 

234.2862 

8 

- 0.3484 

34.4243 

1.1092 

237.7177 

9 

- 0.3084 

34.4338 

1.8727 

239.5300 

10 

- 0.2530 

34.4521 

2.7462 

240.8740 

11 

- 0.1792 

34.4767 

2.4911 

245.2952 

12 

- 0.1191 

34.4967 

3.2346 

249.0272 

13 

- 0.0747 

34.5121 

2.8019 

252.1884 

14 

- 0.0334 

34.5261 

2.2499 

254.8196 

15 

0.0062 

34.5380 

2.3186 

257.1105 

16 

0.0449 

34.5506 

2.3997 

259.7826 

17 

0.0866 

34.5641 

3.5397 

262.5461 

18 

0.1358 

34.5825 

2.3701 

265.8477 

19 

0.1705 

34.5953 

2.8645 

267.9215 


Table 4. One dimensional data set consisting of average 
values of potential temperature, salinity, and layer 
thickness for the entire basin. 


The inputs into the one-dimensional inverse model are 
average values of potential temperature, salinity and layer 
thickness. The total set of 3452 profiles is used to obtain 
averages. These data are shown in Table 4. It is clear 


that temperature and salinity increase with depth, while the 

layer thickness ranges from just over 1 meter to just over 3 
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meters. 


The average thickness is approximately 2.33 meters. 


The average change in temperature between layers is 
approximately 0.046°C, and the average change in salinity 
between layers is approximately 0.16 PSU. 

Model parameters (w, K T , K s ) were calculated for each 
interface above and below layers 2 thru 18. Layers 1 and 19 
are used for input on the boundaries only. Thus 18 values 
of the unknowns are calculated. Vertical velocity (w) is 
found to be constant, which is consistent with the 
integrated form of the continuity equation. As expected, 
heat Flux (F H ) is upward (positive), and K T is negative. 


Interface 

Fh(W/m A 2) 

Kt(10 A -5 m2s A -1) 

w(10 A -7 m/s) 

Rrho 

1 

0.84 

- 0.85 

-5 

6.23 

2 

0.88 

- 0.98 

-5 

5.97 

3 

0.95 

- 1.03 

-5 

5.29 

4 

0.84 

- 1.06 

-5 

5.30 

5 

0.96 

- 1.55 

-5 

5.54 

6 

1.51 

- 1.93 

-5 

5.23 

7 

0.68 

- 2.90 

-5 

7.57 

8 

1.18 

- 1.35 

-5 

3.50 

9 

2.62 

- 1.61 

-5 

4.82 

10 

0.99 

- 1.50 

-5 

4.80 

Mean 

1.56 

-2.31 

-5.00 

5.04 


Table 5. One-dimensional Model 1 results. 


Table 5 is the results from the first one dimensional 
model run. The most important output for this model is the 
Heat Flux (F h ) . The mean heat flux in the upward direction 
is 1.56 W/rrh. This is a critical value to determining 
whether diffusive convection contributes to ice melt. 

The second result is the calculated value of the 
amplitude of C (R p ) . This value is 0.0032 in the Kelley 
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(1990) formulation of the 4/3 flux law (8) . Recall that 
this formulation was derived using laboratory experiments 
that may not represent what is observed in nature. We 
attempt to calibrate this formulation to best represent the 
Arctic data. 
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gure 15. Comparison of the two laboratory representations 
of Heat Flux and the calculated value of Heat Flux from 
the one-dimensional Model 1. 


Figure 15 illustrates that the heat flux (F h ) 
calculated in one-dimensional Model 1 is roughly an order of 
magnitude larger for the Arctic data. 


F h =~K, 


AT 
AZ 


The heat flux was calculated using (36) 

represents the temperature change between the 
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(36) 

Where AT 
layer above 








and layer below the interface, and AZ represents the 
distance between the centers of the two layers. 


Interface 

Fh 

(Model) 

Fh 

(Kelley) 

1 

0.84 

0.13 

2 

0.88 

0.14 

3 

0.95 

0.18 

4 

0.84 

0.20 

5 

0.96 

0.14 

6 

1.51 

0.13 

7 

0.68 

0.04 

8 

1.18 

0.22 

9 

2.62 

0.22 

10 

0.99 

0.33 

Mean 

1.56 

0.18 


Table 6. Comparison of heat flux calculated using the one 
dimensional Model 1 and using Kelley (1990) formulation 

of the 4/3 flux law. 


In this formulation of the heat flux, the magnitude of 
the flux is not only dependent on the change in temperature 
between the two layers, but also on the size of the two 
layers. This proposition is very different from the 4/3 
flux law, which states that the heat flux is independent of 
depth or layer thickness. Comparison with the earlier 
calculations (Kelley, 1990) reveals an order of magnitude 
difference in heat fluxes seen in Table 6. 
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x 10 5 Analysis of fit "Kelley (1990) Fit" for dataset "KtTLS vs. density ratio" 



Figure 16. Fit of C(R p ) as proposed by Kelley (1990) to the 

one dimensional Model 1 results. 

Figure 16 represents an attempt to calibrate the Kelley 
(1990) model to fit the fluxes obtained with the inverse 
model. This resulted in a calculation of the amplitude of 
C (R p ) based on the output of the curve fitting routine. The 
result is a value of the amplitude of C (R p ) = 0.0479, an 
order of magnitude larger than the amplitude of C(R p ) in the 
Kelley (1990) formulation. The upper and lower bounds, at 
the 95% confidence level, are 0.0603 and 0.0354 
respectively. The associated RMSE is 0.0119, and the sum of 
the squared errors is 2.4484e-005. 

The one-dimensional case suggests that the 4/3 flux 
laws as stated by Kelley (1990) and Marmorino and Caldwell 
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(1976) are not applicable to the Arctic diffusive staircases 
in their original form. Instead an order of magnitude 
adjustment must be made to the equations that correspond to 
a new value of C (R p ) . This is similar to the idea of a 
transfer function proposed by Wilson (2007). Additionally, 
the curve is not a perfect fit to the data. Only five of 
the data points lie within the 95% confidence interval. 
There is a much steeper trend associated with this data that 
suggests an even stronger dependence on density ratio. The 
three dimensional model is more inclusive, and shows a trend 
more consistent with the laboratory results. 

B. THREE-DIMENSIONAL MODEL 

1. Depth Dependent Discretization 

Model 1 utilizes the full three-dimensional 
discretization of the equations described by Lee and Veronis 
(1990). The three dimensional data set described in Chapter 
II is used as input to the model. The data set has 425 
individual grid cells for which the model parameters will be 
calculated, resulting in a total of 1275 equations and 1073 
unknowns. The critical values of parameters (w, K s , K T ) 
remain the same as the one dimensional dataset in that they 
are assumed to be constant at each interface. The model 
calculated u and v components of velocity as well. 
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Interface 

Fh(W/m A 2) 

Kt(10 A -5 m2s A -1) 

w(10 A -7 m/s) 

Rrho 

1 

2.15 

3.02 

5.00 

6.36 

2 

2.29 

2.95 

-5.88 

6.07 

3 

2.41 

2.76 

-13.58 

5.29 

4 

1.54 

2.17 

-35.91 

5.50 

5 

1.73 

2.69 

-18.59 

5.66 

6 

1.70 

3.17 

-23.94 

5.44 

7 

3.33 

2.92 

9.82 

5.52 

8 

1.36 

1.88 

-1.67 

4.34 

9 

2.27 

3.01 

29.37 

4.33 

10 

1.91 

2.88 

-6.84 

4.75 

Mean 

2.07 

2.74 

-6.22 

5.32 


Table 7. Three dimensional Model 1 results. 


Table 7 is the results of three-dimensional model 1. 
The heat flux (F h ) is consistent with that of the one 
dimensional model for interfaces 1 through 10. Below 
interface 10 the heat fluxes begin to rise to questionable 
levels. The mean heat flux -2.07 W/m 2 is for the first 10 
interfaces. The assumption of uniform velocity at the 
interface resulted in vastly different values of vertical 
velocity throughout the data. 
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Figure 17. Comparison of the two laboratory representations 
of Heat Flux and the calculated value of Heat Flux from 
the three-dimensional Model 1. 


As in the one dimensional model it is important to show 
how these results relate to the 4/3 flux laws. Figure 17 is 
a plot of the calculated values for heat flux and those 
using the laboratory derived formulations. Like the one 
dimensional model calculations, the three-dimensional model 
values of heat flux are about an order of magnitude larger 
than those calculated using the laboratory formulations of 
the 4/3 flux law. 
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Interface 

Fh 

(Model) 

Fh 

(Kelley) 

1 

2.15 

0.13 

2 

2.29 

0.13 

3 

2.41 

0.15 

4 

1.54 

0.21 

5 

1.73 

0.17 

6 

1.70 

0.10 

7 

3.33 

0.07 

8 

1.36 

0.16 

9 

2.27 

0.26 

10 

1.91 

0.36 

Mean 

2.07 

0.19 


Table 8. Comparison of heat flux calculated using the 
three-dimensional Model 1 and the Kelley (1990) 
formulation of the 4/3 flux law. 

The heat fluxes associated with these values of 
diffusivity are naturally an order of magnitude larger than 
those calculated using the Kelley (1990) formulation. Table 
8 is a comparison of heat flux calculations. The heat flux 
from the model is calculated using (36) . The mean value 
shown for the model calculation is the mean for just the 
first 10 interfaces. The mean for the Kelley formulation is 
that of all 19 interfaces. The difference between the two 
is a full order of magnitude. 
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1.5 


x 10 4 Analysis of fit "Kelley (1990)" for dataset "KtTLS vs. density ratio" 



Figure 18. Fit of C(R p ) as proposed by Kelley (1990) to the 

three dimensional model 1 results. 

Figure 18 is the attempt to calibrate the 4/3 flux law. 
The fit curve excludes the interfaces 14-18 since the 
results for those interfaces were extremely far from the 
rest of the data. The results of the curve fitting were a 
value of the amplitude of C (R p ) = 0.0635. The upper and 
lower bounds of the 95% confidence interval are 0.0805 and 
0.0464 respectively. Nearly all of the data points fall 
within these bounds, which shows that the general trend in 
the data is similar to that of Kelley (1990), yet the model 
results are an order of magnitude larger. 


49 

















2. 4/3 Flux Law Discretization with C(R P ) Unknown 

In order to compare the results of the previous section 
to the 4/3 flux law more directly, the model discretization 
had to be changed. The new model discretization corresponds 
to three-dimensional Model 2 as discussed in Chapter 
III,Section D. 
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The equations remain the same except for the last two 
terms. The coefficients of diffusivity are changed to 
reflect the (AT) 473 relationship described by Turner (1973). 
This relationship is normalized by the mean layer thickness 
and mean temperature gradient in order to maintain the same 
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dimensionality as the original discretization. The 
continuity equation remained the same. The data set, number 
of equations, and number of unknowns are the same as the 
previous section. 

2 fc(h)( A7 T 

F = — - 

i T x 

2h(ATf 

2(c(R p )(&sf 
F s = ~ _, 

2h[ASf 

The equations used to calculate the fluxes are 
different because of the changes to the discretization. 
Temperature (39) and salinity (40) fluxes are now a function 
of (AT) 4/3 , and the fluxes must be normalized by the mean 
layer thickness and mean temperature gradient as well. 


(39) 


(40) 


Interface 

Fh(W/m A 2) 

Kt(10 A -5 m2s A -1) 

w(10 A -7 m/s) 

Rrho 

1 

0.69 

0.78 

5.00 

6.36 

2 

0.79 

0.98 

-12.14 

6.07 

3 

0.86 

1.04 

3.74 

5.29 

4 

0.73 

0.60 

8.01 

5.50 

5 

0.72 

0.73 

-0.90 

5.66 

6 

0.53 

0.96 

13.52 

5.44 

7 

0.45 

1.19 

-1.68 

5.52 

8 

0.56 

0.83 

3.99 

4.34 

9 

0.82 

0.74 

-8.36 

4.33 

10 

0.85 

0.50 

2.47 

4.75 

Mean 

-9.96 

-0.83 

1.37 

5.32 


Table 9. Three-dimensional Model 2 results. 


The model results are presented in Table 9. We again 
observe questionable heat fluxes below layer 10, which we 
attribute to uncertainties of clearly identifying layers 11- 
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15. Those results are not shown. Overall the heat fluxes 
are less than 1 W/rrr. The mean heat flux in the upper 10 
layers is 0.70 W/rrr. 


Interface 

Fh 

(Model) 

Fh 

(Kelley) 

1 

0.69 

0.13 

2 

0.79 

0.13 

3 

0.86 

0.15 

4 

0.73 

0.21 

5 

0.72 

0.17 

6 

0.53 

0.10 

7 

0.45 

0.07 

8 

0.56 

0.16 

9 

0.82 

0.26 

10 

0.85 

0.36 

Mean 

0.70 

0.19 


Table 10. Comparison of heat flux calculated using the 

three-dimensional model 2and the Kelley (1990) 
formulation of the 4/3 flux law. 


Table 10 shows that the fluxes 
calculation are a factor of 3-5 larger 
(1990) calculation. While they are not 
magnitude different as in the previous 
still a significant difference. 


in this inverse 
than the Kelley 
a full order of 
section, there is 
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Figure 19. Comparison of the two laboratory representations 
of C(R p ) and the calculated value of C(R p ) from the 
three-dimensional Model 2. 


Figure 19 is the C(R p ) plot. The general pattern of 
the diffusivities calculated by the inverse model is the 
same. However, the difference between the model and the 
Kelley (1990) and Marmorino and Caldwell (1976) calculations 
is only a factor of 3-5. The five outliers for lower values 
of density ratio correspond to the lower five layers, which 
show anomalous results. 
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Analysis of fit "fit 1" for dataset "KtTLS vs. densityratio" 
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gure 20. Fit of C(R p ) as proposed by Kelley (1990) to the 

three-dimensional Model 2 results. 

Removing the lower five layers from the curve fitting 
routine allows for a much better estimation of the amplitude 
of C (R p ) . Figure 20 is the result of the curve fitting 
routine. The result is a value of the amplitude of C(R p ) = 
0.0185. The upper and lower bounds are equal to 0.0232 and 
0.0139 respectively. This represents a significantly 
reduced confidence interval compared to Model 1. 

3. C(R P ) Calculation using Kelley (1990) 4/3 

Flux Law 

Now we attempt to directly extract the amplitude of the 

coefficient C(R p ) from the inverse calculation. The entire 

Kelley (1990) formulation of the 4/3 flux law is used as the 
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coefficient in the model. Thus, the result is an exact 
value for the amplitude of C (R p ) as model output. We assume 
that the pattern of C (R p ) is correctly captured by the 
laboratory experiments, but the amplitude (B) requires 
calibration. 
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We also assume that B is constant and that it is not 
dependent on the density ratio, so in the model it is 
treated as a single unknown. This reduces the number of 
unknowns to 1038 with the number of equations remaining 
1275. The last two terms are again normalized by the mean 
layer thickness and mean temperature gradient in order to 
maintain the same dimensionality as the original 
discretization. The continuity equation remained the same. 
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Interface 

Fh(W/m A 2) 

w(10 A -7 m/s) 

Rrho 

1 

0.49 

5.00 

6.36 

2 

0.47 

-2.53 

6.07 

3 

0.55 

-5.07 

5.29 

4 

0.79 

11.88 

5.50 

5 

0.63 

-11.01 

5.66 

6 

0.37 

23.70 

5.44 

7 

0.25 

-43.40 

5.52 

8 

0.58 

91.44 

4.34 

9 

0.96 

-108.90 

4.33 

10 

1.34 

82.35 

4.75 

Mean 

0.64 

4.35 

5.32 


Table 11. Three-dimensional Model 3 results. 


The results are shown in Table 11. The heat fluxes are 
generally below 1 W/rrh, with the average heat flux being 
0.70 W/m 2 . The difference is still significant, and is 

comparable to the three dimensional model 2 results. The 
vertical velocities are also anomalously high below layer 
11 . 


Interface 

Fh 

(Model) 

Fh 

(Kelley) 

1 

0.49 

0.13 

2 

0.47 

0.13 

3 

0.55 

0.15 

4 

0.79 

0.21 

5 

0.63 

0.17 

6 

0.37 

0.10 

7 

0.25 

0.07 

8 

0.58 

0.16 

9 

0.96 

0.26 

10 

1.34 

0.36 

Mean 

0.70 

0.19 


Table 12. Comparison of heat flux calculated using the 

three-dimensional Model 3 and the Kelley (1990) 
formulation of the 4/3 flux law. 
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Table 12 is the comparison of Model 3 results and the 
fluxes calculated using the Kelley (1990) formulation of the 
4/3 flux law. The difference between the two is exactly a 
factor of 3.717 throughout the data. This makes sense in 
that all we changed is the amplitude of the calculation, and 
thus all we changed is the multiplying factor determining 
the fluxes. 
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Figure 21. Comparison of the two laboratory representations 
of C (R p ) and the calculated value of C (R p ) from the 
three-dimensional Model 3. 


Figure 21 is the plot of density ratio versus 
diffusivity. The trend in the calculated values from the 
inverse model is the same as the laboratory results. The 
difference in magnitude is approximately a factor of 3. 
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There are no anomalous results since the whole 
parameterization from Kelley (1990) was used. 

The data fitting routine in this case is not necessary, 
since the output of the inverse calculation is a direct 
calculation of the coefficient C. 
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V. DISCUSSION AND CONCLUSIONS 


A. DISCUSSION 

This work attempts to infer the vertical heat fluxes 
through diffusive staircases in the Beaufort Sea using 
techniques of inverse modeling. Table 13 summarizes the 
heat flux results for the four models presented in Chapter 
IV. The values of the fluxes are significantly higher than 
laboratory extrapolation (Kelley, 1990) would suggest. The 
mean fluxes shown are for the upper 10 layers in all of the 
models. The mean of those four is approximately 1.25 W/rrr. 
This suggests that the upward heat flux due to diffusive 
convection is sufficient enough to contribute to sea ice 
melt in the southern Beaufort Sea. 


Interface 

Fh (ID Model) 

Fh (3D Model 1) 

Fh (3D Model 2) 

Fh (3D Model 3) 

Fh (Kelley) 

1 

0.84 

2.15 

0.69 

0.49 

0.13 

2 

0.88 

2.29 

0.79 

0.47 

0.13 

3 

0.95 

2.41 

0.86 

0.55 

0.15 

4 

0.84 

1.54 

0.73 

0.79 

0.21 

5 

0.96 

1.73 

0.72 

0.63 

0.17 

6 

1.51 

1.70 

0.53 

0.37 

0.10 

7 

0.68 

3.33 

0.45 

0.25 

0.07 

8 

1.18 

1.36 

0.56 

0.58 

0.16 

9 

2.62 

2.27 

0.82 

0.96 

0.26 

10 

0.99 

1.91 

0.85 

1.34 

0.36 

Mean 

1.56 

2.07 

0.70 

0.64 

0.19 


Table 13. Comparison of heat fluxes calculated from all 

models. 


The most glaring discrepancies in the calculations are 
the abnormally large flux values below layer 10 in three- 
dimensional models 1 and 2. One reason for these large flux 
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values is the dependence on layer thickness in the model 
discretization. 


Layer 

Thickness (m) 

1 

1.78 

2 

1.59 

3 

1.54 

4 

1.80 

5 

2.12 

6 

2.24 

7 

1.90 

8 

1.14 

9 

1.56 

10 

2.38 

11 

2.45 

12 

3.22 

13 

2.27 

14 

2.67 

15 

1.99 

16 

3.36 

17 

1.72 

18 

2.91 

19 

1.53 


Table 14. Mean layer thicknesses. 


Table 14 is the average layer thicknesses for each 
layer. The deeper the layers the larger the thickness 
values become. There is a significant increase in layer 
thickness below layer 10. Layer thickness is in the 
denominator in the vertical terms. This increase in 
thickness effectively reduces the value of the coefficient 
in the discretized model and requires an increase in the 
calculated unknown. A discretization scheme that uses an 
average layer thickness may reduce this erroneous result. 

The one-dimensional model performs well, and is less 
variable than three-dimensional Models 1 and 2, and its 
variability is comparable to three-dimensional Model 3. 
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Neglecting horizontal velocity likely dampened the effect of 
the layer thickness dependence in the flux calculations. 
Thus, there were no erroneous results in the deeper layers. 


u r = 


8 f dp 


Pof t 


I%* 


v„ =-■ 


8 f d P 


Pof t, dx 




(43) 

(44) 


The vertically integrated thermal wind equations state 
that horizontal velocities are balanced by horizontal 
pressure gradients. This implies that the direction and 
magnitude of the horizontal velocity is strictly controlled 
by the thickness differences between adjacent cells in the 
data. The direction of the flow is calculated to be in the 
direction of the thinner cell, and the magnitude is a result 
of the relative differences in thickness between the two 
cells. 


Layer thickness is in the numerator in the horizontal 
terms, possibly lowering the value of horizontal velocity. 
It is likely that the combination of the sensitivity of the 
vertical heat flux calculations and the sensitivity of the 
horizontal velocity terms to the layer thickness contributed 
to the anomalies in the deeper layers in the two models were 
both were factors in the calculation. 
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Model 

Amplitude of C(Rp) 

1D Model 

0.0479 

3D Model 1 

0.0635 

3D Model 2 

0.0185 

3D Model 3 

0.0119 

Kelley 

0.0032 

Marmorino & Caldwell 

0.0086 


Table 15. Comparison of amplitudes of the 4/3 flux law 

coefficient C (R p ) . 

Table 15 compares the results of each models output for 
the coefficient of amplitude of C (R p ) . As one may expect, 
the formulations of the inverse calculation that more 
closely represents the Kelley (1990) formulation of the 4/3 
flux law (3D Model 2 and 3D Model 3) yielded results that 
were closer to Kelley's value of the amplitude. However, 
these results are still a factor of 4-6 larger than Kelley's 
laboratory results. The one-dimensional model and the 
three-dimensional Model 1 both used the layer thickness 
dependent version of the model formulation and yielded 
results more than a full order of magnitude greater. Caro 
(2009) found results similar to Model 2 and 3 via two- 
dimensional direct numerical simulation, and Wilson (2007) 
found results similar to the ID model and 3D Model 1 
analyzing the ITP data. 

The accuracy of the calculations is not directly 
addressed. There are several questions associated with the 
accuracy of the calculations that need to be answered with 
further study. One very important measure of how well the 
TLS algorithm performs is the condition number. The 
condition number measures the sensitivity of the solution of 
the system of equations to errors in the data. A small (~1) 
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value of the condition number indicates a well conditioned 
system, and thus a system that is not too prone to solution 
error. The larger the condition number, the more sensitive 
the solution is to errors in the system. The condition 
number of the data matrix A in the one dimensional model is 
781 indicating that the system is guite sensitive to data 
errors. The addition of more equations and more unknowns to 
the problem only increases the sensitivity. The condition 
numbers for three-dimensional Models 1-3 are on the order of 
10 8 . 

B. CONCLUSIONS 

The inverse modeling technique was successful at 
calculating the heat flux in the southern Beaufort Sea. 
There are sensitivities associated with the models due to 
the differences in layer thickness. In the lower layers, 
where thickness is greater, there is greater potential for 
large flux values and small horizontal velocities due to the 
larger values of layer thickness. Ultimately, the upper 10 
layers showed that the mean heat fluxes were 1.25 W/m 2 , and 
therefore likely to contribute to sea ice melt. 

The data suggest that the amplitude of the exponential 
form of C (R p ) is likely within the range of a factor of 4 to 
nearly an order of magnitude larger than laboratory results 
(Kelley, 1990; Marmorino & Caldwell, 1976) indicate. The 
application of the inverse model has shown that 
extrapolation of laboratory results to the ocean is not 
perfect, and in this case, not representative of observed 
ocean conditions. Laboratory experiments provide a 
foundation for which observational scientists can pose 
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hypotheses, but ultimately, the true test of a laboratory 
based theory is how well the observed data corroborates the 
results. 

C. RECOMMENDATIONS 

Several questions remain unanswered. First, is there a 
way to improve the condition of the system of equations? 
Second, what is the spatial variability of the heat fluxes 
throuqhout the Beaufort Sea? Finally, can a well desiqned 
experiment, where data is collected systematically in a 
qrid, improve the results? 

It would be beneficial to add a set of equations to 
constrain velocities to known physical principles such as 
Thermal Wind or Conservation of Vorticity. Lee and Veronis 
(1990) did do this, and their results were more consistent 
than estimates based entirely on the advection-diffusion 
equations. 

The spatial variability of the model must be explored. 
This study only used one location in the southern Beaufort 
Sea. The model needs to be run in several locations to qain 
a perspective on the changes in fluxes throughout the 
region. 

The possibility of Inverse Modeling using an 
unstructured grid should be explored. Unstructured modeling 
is a rapidly developing field. This would be very useful in 
oceanography because of the nature of data collection in the 
WHOI ITP program. The use of an unstructured grid would 
allow the direct implementation of these modeling techniques 
without interpolation, thus, reducing one possible source of 
error. 
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In-situ measurements taken within several days would 
likely improve the model. We recommend that in-situ 
measurements need to be collected on a regularly spaced grid 
(similar to C-SALT experiment in Tropics). A clear snapshot 
of the Beaufort Sea at one specific time would significantly 
reduce error in the model. The dynamic nature of the Sea 
creates a source of error naturally. Modeling of this 
nature utilizes static forms of the equations and neglects 
any temporal variation. Thus, a timely set of in-situ 
measurements would likely improve the results. 
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